Covariate descriptive maps

plot_fort_maps2 <- function(data, fill_str, legend_title, plot_title, size=0.1){
  require(ggplot2)
  require(ggsn)
  require(broom)
  ggplot() +                                               # initialize ggplot object
    geom_polygon(                                          # make a polygon
      data = data,                                    # data frame
      aes_string(x = "long", y = "lat", group = "group",                # coordinates, and group them by polygons
                 fill = fill_str),
      size=size, color="black") +                # variable to use for filling
    scale_fill_brewer(name=legend_title, palette = "RdYlBu", direction = -1,
                      drop = FALSE) + # fill with brewer colors   # add title
    theme(line = element_blank(),
          axis.text=element_blank(),      # .. tickmarks..
          axis.title=element_blank(),
          legend.position="none",
          plot.margin = unit(c(0, 0, 0, 0), "cm"),
          #legend.text=element_text(size=15),
          #legend.title=element_text(size=17), # .. axis labels..
          panel.background = element_blank(), plot.title = element_text(size=5)) + ggtitle(plot_title)
}
shape_to_ggplot <- function(shape){
  require(broom)
  gg_data <- tidy(shape)
  data <- slot(shape, "data")
  shape[["polyID"]] <- sapply(slot(shape, "polygons"), function(x) slot(x, "ID"))
  gg_data <- merge(gg_data, shape, by.x="id", by.y="polyID")
  return(gg_data)
}
grid_

Error: object 'grid_' not found

library(readr)
library(dplyr)
library(stringi)
library(rgdal)
CT_shp <- readOGR("~/Documents/Harvard - SM80/Thesis/Census_data/CE_setores_censitarios_shape/", "Corrected_CTs", use_iconv = TRUE, encoding = "latin1")

OGR data source with driver: ESRI Shapefile 
Source: "/Users/Sudipta/Documents/Harvard - SM80/Thesis/Census_data/CE_setores_censitarios_shape", layer: "Corrected_CTs"
with 3044 features
It has 13 fields

library(spdep)

package ‘spdep’ was built under R version 3.4.4Loading required package: Matrix
Loading required package: spData
package ‘spData’ was built under R version 3.4.4To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source'))`

library(dplyr)
CT_shp@data$CD_GEOCODI <- as.character(CT_shp@data$CD_GEOCODI)
CT_shp@data <- arrange(CT_shp@data, CD_GEOCODI)

package ‘bindrcpp’ was built under R version 3.4.4

CT_wm <- poly2nb(CT_shp, snap=0.001)

library(INLA)

This is INLA_17.06.20 built 2017-06-20 03:44:36 UTC.
See www.r-inla.org/contact-us for how to get help.

nb2INLA("CT_graph.adj", CT_wm)
CT_adj <- "CT_graph.adj"

Census_covar <- read_csv("~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit/Census_data/census_covariates_Fortaleza_01_03_19.csv")

Parsed with column specification:
cols(
  .default = col_integer(),
  Cod_setor = col_double(),
  TIPO = col_character(),
  CD_GEOCODB = col_double(),
  NM_BAIRRO = col_character(),
  Lit_rate = col_double(),
  Mean_HH_inc = col_double(),
  Perc_Wtr_Spply_Ntwrk = col_double(),
  Perc_Sewerage_Conn = col_double(),
  Perc_Garbage_Col_Serv = col_double(),
  Perc_PPH_Elec = col_double(),
  Perc_branca = col_double(),
  Perc_preta = col_double(),
  Perc_amarela = col_double(),
  Perc_parda = col_double(),
  Perc_indigena = col_double(),
  ICE = col_double()
)
See spec(...) for full column specifications.

Census_covar[is.na(Census_covar$Mean_HH_inc),]$Mean_HH_inc <- 1
Census_covar[Census_covar$Mean_HH_inc==0,]$Mean_HH_inc <- 1
Census_covar$log_Mean_HH_inc <- log(Census_covar$Mean_HH_inc)
Census_covar$Cod_setor <- as.character(Census_covar$Cod_setor)
Census_covar$LII <- NA
for (i in 1:3044){
  CT_Code <- CT_shp@data$CD_GEOCODI[i]
  neighbors_codes <- CT_shp@data$CD_GEOCODI[CT_wm[[i]]]
  total_nb_inc <- sum(Census_covar[(Census_covar$Cod_setor %in% neighbors_codes),]$Total_HH_inc_in_CT, na.rm = TRUE)
  total_nb_hh <- sum(Census_covar[(Census_covar$Cod_setor %in% neighbors_codes),]$Total_no_HH_in_CT, na.rm = TRUE)
  Mean_nb_inc <- total_nb_inc / total_nb_hh
  LII <- Mean_nb_inc / Census_covar[(Census_covar$Cod_setor)==CT_Code,]$Mean_HH_inc
  Census_covar[(Census_covar$Cod_setor==CT_Code),]$LII <- LII
}
Census_covar$LII[Census_covar$LII==Inf] <- NA 
Census_covar$log_LII <- log(Census_covar$LII)

CT_shp_gg <- shape_to_ggplot(CT_shp)

Loading required package: broom
package ‘broom’ was built under R version 3.4.4Regions defined for each Polygons

breaks_inc <- quantile(Census_covar$log_Mean_HH_inc, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
breaks_lii <- quantile(Census_covar$LII, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
breaks_lit <- quantile(Census_covar$Lit_rate, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
breaks_log_lii <- quantile(Census_covar$log_LII, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
breaks_brnc <- quantile(Census_covar$Perc_branca, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
breaks_swrge <- quantile(Census_covar$Perc_Sewerage_Conn, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
breaks_ice <- quantile(Census_covar$ICE, probs=c(seq(0,1,0.15),1), na.rm = TRUE)
CT_shp_gg_covar <- CT_shp_gg[,-1] %>% left_join(Census_covar, by=c("CD_GEOCODI"="Cod_setor")) %>% mutate(cut_inc=cut(log_Mean_HH_inc, breaks=breaks_inc, include.lowest=TRUE), cut_lii=cut(LII, breaks=breaks_lii, include.lowest=TRUE), cut_brnc=cut(Perc_branca, breaks=breaks_brnc, include.lowest=TRUE), cut_swrge=cut(Perc_Sewerage_Conn, breaks=breaks_swrge, include.lowest=TRUE), cut_ice=cut(ICE, breaks=breaks_ice, include.lowest=TRUE), cut_log_lii=cut(log_LII, breaks=breaks_log_lii, include.lowest=TRUE), cut_lit=cut(Lit_rate, breaks=breaks_lit, include.lowest=TRUE))

plot_fort_maps2(CT_shp_gg_covar, "cut_inc","Log Mean HH Income", "Log Mean HH Income") + theme(plot.title = element_text(size=10), legend.position = "right") + scale_fill_brewer(name="Log Mean HH Income", palette = "RdYlBu", direction=1, drop = FALSE)

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing
scale.

plot_fort_maps2(CT_shp_gg_covar, "cut_log_lii","Log LII", "Log LII") + theme(plot.title = element_text(size=10), legend.position = "right") 

plot_fort_maps2(CT_shp_gg_covar, "cut_lii","LII", "LII") + theme(plot.title = element_text(size=10), legend.position = "right") 

plot_fort_maps2(CT_shp_gg_covar, "cut_lit","Literacy Rate", "Literacy Rate") + theme(plot.title = element_text(size=10), legend.position = "right") + scale_fill_brewer(name="Literacy Rate", palette = "RdYlBu", direction=1, drop = FALSE)

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing
scale.

plot_fort_maps2(CT_shp_gg_covar, "cut_brnc","Perc branca", "Perc branca") + theme(plot.title = element_text(size=10), legend.position = "right") + scale_fill_brewer(name="Perc branca", palette = "RdYlBu", direction=1, drop = FALSE)

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing
scale.

plot_fort_maps2(CT_shp_gg_covar, "cut_swrge","Perc Sewerage", "Perc Sewerage") + theme(plot.title = element_text(size=10), legend.position = "right") + scale_fill_brewer(name="Perc Sewerage", palette = "RdYlBu", direction=1, drop = FALSE)

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing
scale.

plot_fort_maps2(CT_shp_gg_covar, "cut_ice","ICE Index", "ICE Index") + theme(plot.title = element_text(size=10), legend.position = "right") + scale_fill_brewer(name="ICE Index", palette = "RdYlBu", direction=1, drop = FALSE)

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing
scale.